]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C
new task MK
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawCorrelationFunctionPsi.C
index fe4a157597c6bb663086abfec42453c1559f83e6..fd5f0954934402c4b0bd6b920908caf9882939eb 100644 (file)
@@ -1,9 +1,10 @@
 const Int_t numberOfCentralityBins = 12;
-TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
+//TString centralityArray[numberOfCentralityBins] = {"0-4","4-5","6-14","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
+TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
 
 const Int_t gRebin = 1;
 
-void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root", 
+void drawCorrelationFunctionPsi(const char* filename = "AnalysisResultsPsi_train_06feb.root", 
                                Int_t gCentrality = 1,
                                Int_t gBit = -1,
                                const char* gCentralityEstimator = 0x0,
@@ -11,11 +12,14 @@ void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root",
                                Bool_t kShowMixed = kTRUE, 
                                Double_t psiMin = -0.5, 
                                Double_t psiMax = 3.5,
+                               Double_t vertexZMin = -10.,
+                               Double_t vertexZMax = 10.,
                                Double_t ptTriggerMin = -1.,
                                Double_t ptTriggerMax = -1.,
                                Double_t ptAssociatedMin = -1.,
                                Double_t ptAssociatedMax = -1.,
-                               Bool_t normToTrig = kFALSE,
+                               Bool_t normToTrig = kTRUE,
+                               Bool_t kUseVzBinning = kTRUE,
                                Int_t rebinEta = 1,
                                Int_t rebinPhi = 1,
                                TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
@@ -23,11 +27,20 @@ void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root",
   //Macro that draws the correlation functions from the balance function
   //analysis vs the reaction plane
   //Author: Panos.Christakoglou@nikhef.nl
-  gROOT->LoadMacro("~/SetPlotStyle.C");
-  SetPlotStyle();
+  //gROOT->LoadMacro("~/SetPlotStyle.C");
+  //SetPlotStyle();
   gStyle->SetPalette(1,0);
 
   //Load the PWG2ebye library
+  gSystem->Load("libCore.so");        
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libTree.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so");
+  
   gSystem->Load("libANALYSIS.so");
   gSystem->Load("libANALYSISalice.so");
   gSystem->Load("libEventMixing.so");
@@ -51,6 +64,7 @@ void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root",
       Printf("The file %s is not found.",filename);
     }
     
+    
     TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
     if(!dir) {   
       Printf("The TDirectoryFile is not found.",filename);
@@ -76,8 +90,8 @@ void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root",
   }
   else 
     draw(list,listShuffled,listMixed,listQA,
-        gCentralityEstimator,gCentrality,psiMin,psiMax,
-        ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig,rebinEta,rebinPhi,eventClass);
+        gCentralityEstimator,gCentrality,psiMin,psiMax,vertexZMin,vertexZMax,
+        ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig,kUseVzBinning,rebinEta,rebinPhi,eventClass);
 }
 
 //______________________________________________________//
@@ -103,7 +117,7 @@ TList *GetListOfObjects(const char* filename,
     return listBF;
   }
   //dir->ls();
-  
+
   TString listBFName;
   if(kData == 0) {
     //cout<<"no shuffling - no mixing"<<endl;
@@ -240,12 +254,16 @@ TList *GetListOfObjects(const char* filename,
 }
 
 //______________________________________________________//
-void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
+void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, 
+         TList *listQA,
          const char *gCentralityEstimator,
-         Int_t gCentrality, Double_t psiMin, Double_t psiMax,
+         Int_t gCentrality, Double_t psiMin, Double_t psiMax,  
+         Double_t vertexZMin,Double_t vertexZMax,
          Double_t ptTriggerMin, Double_t ptTriggerMax,
          Double_t ptAssociatedMin, Double_t ptAssociatedMax,
-         Bool_t normToTrig, Int_t rebinEta, Int_t rebinPhi,TString eventClass) {
+         Bool_t normToTrig,                            
+         Bool_t kUseVzBinning,
+         Int_t rebinEta, Int_t rebinPhi,TString eventClass) {
   //Draws the correlation functions for every centrality bin
   //(+-), (-+), (++), (--)  
   AliTHn *hP = NULL;
@@ -286,6 +304,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   b->SetHistNnp(hNP);
   b->SetHistNpp(hPP);
   b->SetHistNnn(hNN);
+  if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
 
   //balance function shuffling
   AliTHn *hPShuffled = NULL;
@@ -331,6 +350,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     bShuffled->SetHistNnp(hNPShuffled);
     bShuffled->SetHistNpp(hPPShuffled);
     bShuffled->SetHistNnn(hNNShuffled);
+    if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
   }
 
   //balance function mixing
@@ -377,6 +397,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     bMixed->SetHistNnp(hNPMixed);
     bMixed->SetHistNpp(hPPMixed);
     bMixed->SetHistNnn(hNNMixed);
+    if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
   }
 
 
@@ -429,7 +450,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " % ";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else if(eventClass == "Multiplicity"){
     histoTitle = "(+-) | Multiplicity: ";
@@ -437,23 +458,22 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " tracks";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else{ // "EventPlane" (default)
     histoTitle = "(+-) | Centrality: ";
     histoTitle += centralityArray[gCentrality-1]; 
     histoTitle += "%";
     if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
+      histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})"; 
     else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
+      histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})"; 
     else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
+      histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})"; 
     else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+      histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
-
-  gHistPN[0] = b->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+  gHistPN[0] = b->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
   if(rebinEta > 1 || rebinPhi > 1){
     gHistPN[0]->Rebin2D(rebinEta,rebinPhi);
     gHistPN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -475,19 +495,10 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   //cPN[0]->SaveAs(pngName.Data());
   
   if(listBFShuffled) {
-    histoTitle = "(+-) shuffled | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    histoTitle.ReplaceAll("(+-)","(+-) shuffled"); 
     
-    gHistPN[1] = bShuffled->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistPN[1] = bShuffled->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPN[1]->Rebin2D(rebinEta,rebinPhi);
       gHistPN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -511,20 +522,12 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
 
   if(listBFMixed) {
-    histoTitle = "(+-) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) mixed"); 
+    else                                histoTitle.ReplaceAll("(+-)","(+-) mixed"); 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
       gHistPN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -555,10 +558,13 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     //cPN[2]->SaveAs(pngName.Data());
 
     //Correlation function (+-)
-    gHistPN[3] = dynamic_cast<TH2D *>(gHistPN[0]->Clone());
-    gHistPN[3]->Divide(gHistPN[2]);
+    gHistPN[3] = b->GetCorrelationFunction("PN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
     gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+    //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
     cPN[3] = new TCanvas("cPN3","",0,300,600,500);
     cPN[3]->SetFillColor(10); 
     cPN[3]->SetHighLightColor(10);
@@ -575,17 +581,9 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
   // if no mixing then divide by convoluted histograms
   else if(listQA){
-    histoTitle = "(+-) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) convoluted"); 
+    else                                histoTitle.ReplaceAll("(+-)","(+-) convoluted"); 
     
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
@@ -620,7 +618,11 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     gHistPN[3] = dynamic_cast<TH2D *>(gHistPN[0]->Clone());
     gHistPN[3]->Divide(gHistPN[2]);
     gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+    //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
     cPN[3] = new TCanvas("cPN3","",0,300,600,500);
     cPN[3]->SetFillColor(10); 
     cPN[3]->SetHighLightColor(10);
@@ -643,7 +645,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " % ";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else if(eventClass == "Multiplicity"){
     histoTitle = "(-+) | Multiplicity: ";
@@ -651,23 +653,23 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " tracks";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else{ // "EventPlane" (default)
     histoTitle = "(-+) | Centrality: ";
     histoTitle += centralityArray[gCentrality-1]; 
     histoTitle += "%";
     if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
+      histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})"; 
     else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
+      histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})"; 
     else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
+      histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})"; 
     else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+      histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
 
-  gHistNP[0] = b->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+  gHistNP[0] = b->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
   if(rebinEta > 1 || rebinPhi > 1){
     gHistNP[0]->Rebin2D(rebinEta,rebinPhi);
     gHistNP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -690,19 +692,10 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   //cNP[0]->SaveAs(pngName.Data());
 
   if(listBFShuffled) {
-    histoTitle = "(-+) shuffled | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    histoTitle.ReplaceAll("(-+)","(-+) shuffled"); 
     
-    gHistNP[1] = bShuffled->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistNP[1] = bShuffled->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNP[1]->Rebin2D(rebinEta,rebinPhi);
       gHistNP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -726,20 +719,12 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
 
   if(listBFMixed) {
-    histoTitle = "(-+) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
-
+    
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) mixed"); 
+    else                                histoTitle.ReplaceAll("(-+)","(-+) mixed"); 
+    
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
       gHistNP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
@@ -769,10 +754,13 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     //cNP[2]->SaveAs(pngName.Data());
 
     //Correlation function (-+)
-    gHistNP[3] = dynamic_cast<TH2D *>(gHistNP[0]->Clone());
-    gHistNP[3]->Divide(gHistNP[2]);
+    gHistNP[3] = b->GetCorrelationFunction("NP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
     gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+    //gHistNP[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
     cNP[3] = new TCanvas("cNP3","",100,300,600,500);
     cNP[3]->SetFillColor(10); 
     cNP[3]->SetHighLightColor(10);
@@ -789,17 +777,9 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
   // if no mixing then divide by convoluted histograms
   else if(listQA){
-    histoTitle = "(-+) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) convoluted"); 
+    else                                histoTitle.ReplaceAll("(-+)","(-+) convoluted"); 
 
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
@@ -834,7 +814,11 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     gHistNP[3] = dynamic_cast<TH2D *>(gHistNP[0]->Clone());
     gHistNP[3]->Divide(gHistNP[2]);
     gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+    //gHistNP[3]->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
     cNP[3] = new TCanvas("cNP3","",100,300,600,500);
     cNP[3]->SetFillColor(10); 
     cNP[3]->SetHighLightColor(10);
@@ -858,7 +842,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " % ";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else if(eventClass == "Multiplicity"){
     histoTitle = "(++) | Multiplicity: ";
@@ -866,23 +850,23 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " tracks";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else{ // "EventPlane" (default)
     histoTitle = "(++) | Centrality: ";
     histoTitle += centralityArray[gCentrality-1]; 
     histoTitle += "%";
     if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
+      histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})"; 
     else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
+      histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})"; 
     else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
+      histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})"; 
     else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+      histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
 
-  gHistPP[0] = b->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+  gHistPP[0] = b->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
   if(rebinEta > 1 || rebinPhi > 1){
     gHistPP[0]->Rebin2D(rebinEta,rebinPhi);
     gHistPP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -905,19 +889,10 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   //cPP[0]->SaveAs(pngName.Data());
   
   if(listBFShuffled) {
-    histoTitle = "(++) shuffled | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    histoTitle.ReplaceAll("(++)","(++) shuffled"); 
     
-    gHistPP[1] = bShuffled->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistPP[1] = bShuffled->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPP[1]->Rebin2D(rebinEta,rebinPhi);
       gHistPP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -941,20 +916,12 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
 
   if(listBFMixed) {
-    histoTitle = "(++) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) mixed"); 
+    else                                histoTitle.ReplaceAll("(++)","(++) mixed"); 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
       gHistPP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -984,10 +951,13 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     //cPP[2]->SaveAs(pngName.Data());
 
     //Correlation function (++)
-    gHistPP[3] = dynamic_cast<TH2D *>(gHistPP[0]->Clone());
-    gHistPP[3]->Divide(gHistPP[2]);
+    gHistPP[3] = b->GetCorrelationFunction("PP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
     gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+    // gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
     cPP[3] = new TCanvas("cPP3","",200,300,600,500);
     cPP[3]->SetFillColor(10); 
     cPP[3]->SetHighLightColor(10);
@@ -1004,17 +974,9 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
   // if no mixing then divide by convoluted histograms
   else if(listQA){
-        histoTitle = "(++) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) convoluted"); 
+    else                                histoTitle.ReplaceAll("(++)","(++) convoluted"); 
     
     if(rebinEta > 1 || rebinPhi > 1){
       gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
@@ -1048,7 +1010,11 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     gHistPP[3] = dynamic_cast<TH2D *>(gHistPP[0]->Clone());
     gHistPP[3]->Divide(gHistPP[2]);
     gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+    //gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
     cPP[3] = new TCanvas("cPP3","",200,300,600,500);
     cPP[3]->SetFillColor(10); 
     cPP[3]->SetHighLightColor(10);
@@ -1071,7 +1037,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " % ";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else if(eventClass == "Multiplicity"){
     histoTitle = "(--) | Multiplicity: ";
@@ -1079,23 +1045,23 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     histoTitle += " - ";
     histoTitle += psiMax;
     histoTitle += " tracks";
-    histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   }
   else{ // "EventPlane" (default)
     histoTitle = "(--) | Centrality: ";
     histoTitle += centralityArray[gCentrality-1]; 
     histoTitle += "%";
     if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
+      histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})"; 
     else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
+      histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})"; 
     else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
+      histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})"; 
     else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+      histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})"; 
   } 
 
-  gHistNN[0] = b->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+  gHistNN[0] = b->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
   if(rebinEta > 1 || rebinPhi > 1){
     gHistNN[0]->Rebin2D(rebinEta,rebinPhi);
     gHistNN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -1118,19 +1084,10 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   //cNN[0]->SaveAs(pngName.Data());
 
   if(listBFShuffled) {
-    histoTitle = "(--) shuffled | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    histoTitle.ReplaceAll("(--)","(--) shuffled"); 
     
-    gHistNN[1] = bShuffled->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistNN[1] = bShuffled->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNN[1]->Rebin2D(rebinEta,rebinPhi);
       gHistNN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -1154,20 +1111,12 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
 
   if(listBFMixed) {
-    histoTitle = "(--) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) mixed"); 
+    else                                histoTitle.ReplaceAll("(--)","(--) mixed"); 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
+    gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
       gHistNN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));  
@@ -1197,10 +1146,13 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     //cNN[2]->SaveAs(pngName.Data());
 
     //Correlation function (--)
-    gHistNN[3] = dynamic_cast<TH2D *>(gHistNN[0]->Clone());
-    gHistNN[3]->Divide(gHistNN[2]);
+    gHistNN[3] = b->GetCorrelationFunction("NN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
     gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+    // gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
     cNN[3] = new TCanvas("cNN3","",300,300,600,500);
     cNN[3]->SetFillColor(10); 
     cNN[3]->SetHighLightColor(10);
@@ -1217,17 +1169,9 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   }
   // if no mixing then divide by convoluted histograms
   else if(listQA){
-        histoTitle = "(--) mixed | Centrality: "; 
-    histoTitle += centralityArray[gCentrality-1]; 
-    histoTitle += "%";
-    if((psiMin == -0.5)&&(psiMax == 0.5))
-      histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})"; 
-    else if((psiMin == 0.5)&&(psiMax == 1.5))
-      histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})"; 
-    else if((psiMin == 1.5)&&(psiMax == 2.5))
-      histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})"; 
-    else 
-      histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
+
+    if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) convoluted"); 
+    else                                histoTitle.ReplaceAll("(--)","(--) convoluted"); 
     
     if(rebinEta > 1 || rebinPhi > 1){
       gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
@@ -1261,7 +1205,11 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
     gHistNN[3] = dynamic_cast<TH2D *>(gHistNN[0]->Clone());
     gHistNN[3]->Divide(gHistNN[2]);
     gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
-    gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+    if(normToTrig)
+      gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
+    else
+      gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+    //gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
     cNN[3] = new TCanvas("cNN3","",300,300,600,500);
     cNN[3]->SetFillColor(10); 
     cNN[3]->SetHighLightColor(10);
@@ -1300,7 +1248,13 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
   newFileName += Form("%.1f",ptAssociatedMax); 
+
+  newFileName += "_"; 
+  newFileName += Form("%.1f",psiMin);
+  newFileName += "-";   
+  newFileName += Form("%.1f",psiMax);
   newFileName += ".root";
+
   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
   gHistPN[0]->SetName("gHistPNRaw"); gHistPN[0]->Write();
   gHistNP[0]->SetName("gHistNPRaw"); gHistNP[0]->Write();
@@ -1366,14 +1320,19 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed, TList *listQA,
 }
 
 //____________________________________________________________//
-void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
-                             Int_t gTrainID = 208,                           
+void drawCorrelationFunctions(const char* lhcPeriod = "LHC10h",
+                             const char* gCentralityEstimator = "V0M",
+                             Int_t gBit = 128,
+                             const char* gEventPlaneEstimator = "VZERO",
                              Int_t gCentrality = 1,
                              Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+                             Double_t vertexZMin = -10.,
+                             Double_t vertexZMax = 10.,
                              Double_t ptTriggerMin = -1.,
                              Double_t ptTriggerMax = -1.,
                              Double_t ptAssociatedMin = -1.,
-                             Double_t ptAssociatedMax = -1.) {
+                             Double_t ptAssociatedMax = -1.,
+                             Bool_t kFit = kFALSE) {
   //Macro that draws the charge dependent correlation functions
   //for each centrality bin for the different pT of trigger and 
   //associated particles
@@ -1381,10 +1340,17 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   TGaxis::SetMaxDigits(3);
 
   //Get the input file
-  TString filename = "PbPb/"; filename += lhcPeriod; 
-  filename +="/Train"; filename += gTrainID;
-  filename +="/Centrality"; filename += gCentrality;
-  filename += "/correlationFunction.Centrality";
+  /* TString filename = lhcPeriod; 
+  filename += "/Centrality"; filename += gCentralityEstimator;
+  filename += "_Bit"; filename += gBit;
+  filename += "_"; filename += gEventPlaneEstimator;
+  filename +="/PttFrom";
+  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
+  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
+  filename += Form("%.1f",ptAssociatedMax); */
+
+  TString filename = "correlationFunction.Centrality";
   filename += gCentrality; filename += ".Psi";
   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
@@ -1395,13 +1361,18 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
   filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
   filename += Form("%.1f",ptAssociatedMax); 
+  
+  filename += "_"; 
+  filename += Form("%.1f",psiMin);
+  filename += "-";   
+  filename += Form("%.1f",psiMax);
   filename += ".root";  
 
   //Open the file
   TFile *f = TFile::Open(filename.Data());
   if((!f)||(!f->IsOpen())) {
     Printf("The file %s is not found. Aborting...",filename);
-    return listBF;
+    return;
   }
   //f->ls();
   
@@ -1412,13 +1383,13 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
 
   TString psiLatex;
   if((psiMin == -0.5)&&(psiMax == 0.5))
-    psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}"; 
+    psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}"; 
   else if((psiMin == 0.5)&&(psiMax == 1.5))
-    psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}"; 
+    psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}"; 
   else if((psiMin == 1.5)&&(psiMax == 2.5))
-    psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}"; 
+    psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}"; 
   else 
-    psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}"; 
+    psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}"; 
  
   TString pttLatex = Form("%.1f",ptTriggerMin);
   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
@@ -1439,25 +1410,26 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   //Get the +- correlation function
   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
   gHistPN->SetStats(kFALSE);
-  gHistPN->SetTitle("");
-  gHistPN->GetXaxis()->SetRangeUser(-1.45,1.45);
+  gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+  gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
   gHistPN->GetXaxis()->CenterTitle();
   gHistPN->GetXaxis()->SetTitleOffset(1.2);
   gHistPN->GetYaxis()->CenterTitle();
   gHistPN->GetYaxis()->SetTitleOffset(1.2);
-  gHistPN->GetZaxis()->SetTitleOffset(1.2);
+  gHistPN->GetZaxis()->SetTitleOffset(1.5);
   TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
   cPN->SetFillColor(10); cPN->SetHighLightColor(10);
   cPN->SetLeftMargin(0.15);
   gHistPN->DrawCopy("surf1fb");
+
   gPad->SetTheta(30); // default is 30
   gPad->SetPhi(-60); // default is 30
   gPad->Update();
 
-  latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
-  //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
 
   pngName = "CorrelationFunction.Centrality"; 
   pngName += centralityArray[gCentrality-1]; 
@@ -1473,32 +1445,35 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   pngName += Form("%.1f",ptAssociatedMax); 
   pngName += ".PositiveNegative.png";
   cPN->SaveAs(pngName.Data());
-  fitCorrelationFunctions(gCentrality, psiMin, psiMax,
-                         ptTriggerMin,ptTriggerMax,
-                         ptAssociatedMin, ptAssociatedMax,gHistPN);
+  if(kFit)
+    fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
+                           ptTriggerMin,ptTriggerMax,
+                           ptAssociatedMin, ptAssociatedMax,gHistPN);
+
   //============================================================//
   //Get the -+ correlation function
   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
   gHistNP->SetStats(kFALSE);
-  gHistNP->SetTitle("");
-  gHistNP->GetXaxis()->SetRangeUser(-1.45,1.45);
+  gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+  gHistNP->GetXaxis()->SetRangeUser(-1.4,1);
   gHistNP->GetXaxis()->CenterTitle();
   gHistNP->GetXaxis()->SetTitleOffset(1.2);
   gHistNP->GetYaxis()->CenterTitle();
   gHistNP->GetYaxis()->SetTitleOffset(1.2);
-  gHistNP->GetZaxis()->SetTitleOffset(1.2);
+  gHistNP->GetZaxis()->SetTitleOffset(1.5);
   TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
   cNP->SetFillColor(10); cNP->SetHighLightColor(10);
   cNP->SetLeftMargin(0.15);
   gHistNP->DrawCopy("surf1fb");
+
   gPad->SetTheta(30); // default is 30
   gPad->SetPhi(-60); // default is 30
   gPad->Update();
 
-  latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
-  //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
 
   pngName = "CorrelationFunction.Centrality"; 
   pngName += centralityArray[gCentrality-1]; 
@@ -1515,32 +1490,35 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   pngName += ".NegativePositive.png";
   cNP->SaveAs(pngName.Data());
 
-  fitCorrelationFunctions(gCentrality, psiMin, psiMax,
-                         ptTriggerMin,ptTriggerMax,
-                         ptAssociatedMin, ptAssociatedMax,gHistNP);
+  if(kFit)
+    fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
+                           ptTriggerMin,ptTriggerMax,
+                           ptAssociatedMin, ptAssociatedMax,gHistNP);
+
   //============================================================//
   //Get the ++ correlation function
   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
   gHistPP->SetStats(kFALSE);
-  gHistPP->SetTitle("");
-  gHistPP->GetXaxis()->SetRangeUser(-1.45,1.45);
+  gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+  gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
   gHistPP->GetXaxis()->CenterTitle();
   gHistPP->GetXaxis()->SetTitleOffset(1.2);
   gHistPP->GetYaxis()->CenterTitle();
   gHistPP->GetYaxis()->SetTitleOffset(1.2);
-  gHistPP->GetZaxis()->SetTitleOffset(1.2);
+  gHistPP->GetZaxis()->SetTitleOffset(1.5);
   TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
   cPP->SetFillColor(10); cPP->SetHighLightColor(10);
   cPP->SetLeftMargin(0.15);
   gHistPP->DrawCopy("surf1fb");
+
   gPad->SetTheta(30); // default is 30
   gPad->SetPhi(-60); // default is 30
   gPad->Update();
 
-  latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
-  //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
 
   pngName = "CorrelationFunction.Centrality"; 
   pngName += centralityArray[gCentrality-1]; 
@@ -1557,32 +1535,35 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   pngName += ".PositivePositive.png";
   cPP->SaveAs(pngName.Data());
 
-  fitCorrelationFunctions(gCentrality, psiMin, psiMax,
-                         ptTriggerMin,ptTriggerMax,
-                         ptAssociatedMin, ptAssociatedMax,gHistPP);
+  if(kFit)
+    fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
+                           ptTriggerMin,ptTriggerMax,
+                           ptAssociatedMin, ptAssociatedMax,gHistPP);
+
   //============================================================//
   //Get the -- correlation function
   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
   gHistNN->SetStats(kFALSE);
-  gHistNN->SetTitle("");
-  gHistNN->GetXaxis()->SetRangeUser(-1.45,1.45);
+  gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+  gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
   gHistNN->GetXaxis()->CenterTitle();
   gHistNN->GetXaxis()->SetTitleOffset(1.2);
   gHistNN->GetYaxis()->CenterTitle();
   gHistNN->GetYaxis()->SetTitleOffset(1.2);
-  gHistNN->GetZaxis()->SetTitleOffset(1.2);
+  gHistNN->GetZaxis()->SetTitleOffset(1.5);
   TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
   cNN->SetFillColor(10); cNN->SetHighLightColor(10);
   cNN->SetLeftMargin(0.15);
   gHistNN->DrawCopy("surf1fb");
+
   gPad->SetTheta(30); // default is 30
   gPad->SetPhi(-60); // default is 30
   gPad->Update();
 
-  latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
-  //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
-  latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
 
   pngName = "CorrelationFunction.Centrality"; 
   pngName += centralityArray[gCentrality-1]; 
@@ -1599,150 +1580,559 @@ void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
   pngName += ".NegativeNegative.png";
   cNN->SaveAs(pngName.Data());
 
-  fitCorrelationFunctions(gCentrality, psiMin, psiMax,
-                         ptTriggerMin,ptTriggerMax,
-                         ptAssociatedMin, ptAssociatedMax,gHistNN);
+  if(kFit)
+    fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
+                           ptTriggerMin,ptTriggerMax,
+                           ptAssociatedMin, ptAssociatedMax,gHistNN);
 }
 
-// //____________________________________________________________//
-// void fitCorrelationFunctions(Int_t gCentrality = 1,
-//                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
-//                          Double_t ptTriggerMin = -1.,
-//                          Double_t ptTriggerMax = -1.,
-//                          Double_t ptAssociatedMin = -1.,
-//                          Double_t ptAssociatedMax = -1.,
-//                          TH2D *gHist) {
-
-//   cout<<"FITTING FUNCTION (MW style)"<<endl;
-
-//   //near side peak(s): [1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) + 
-//   //                        TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4])))
-//   //away side peak(s): [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power((y/[9]),2)),[10])) + 
-//   //                        TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power((y/[9]),2)),[10])))
-//   //flow contribution (v1 up to v4): 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))
-
-
-
-//   TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4]))) + [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10]))) + 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.); 
-//   gFitFunction->SetName("gFitFunction");
-
-
-//   //Normalization
-//   gFitFunction->SetParName(0,"N1"); 
-//   //near side peak(s)
-//   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
-//   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(2,0.0);
-//   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
-//   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
-//   gFitFunction->SetParName(5,"Offset_{near side}"); gFitFunction->FixParameter(5,0.0);
-//   gFitFunction->SetParName(6,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(6,0.0);
-
-//   //away side peak(s)
-//   gFitFunction->SetParName(7,"N_{near side}");gFitFunction->FixParameter(7,0.0);
-//   gFitFunction->SetParName(8,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(8,0.0);
-//   gFitFunction->SetParName(9,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(9,0.0);
-//   gFitFunction->SetParName(10,"Exponent_{near side}"); gFitFunction->FixParameter(10,1.0);
-//   gFitFunction->SetParName(11,"Offset_{near side}"); gFitFunction->FixParameter(11,0.0);
-//   gFitFunction->SetParName(12,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(12,0.0);
-
-//   //flow contribution
-//   gFitFunction->SetParName(13,"N_{flow}"); gFitFunction->SetParameter(13,0.2);
-//   gFitFunction->SetParName(14,"V1"); gFitFunction->SetParameter(14,0.005);
-//   gFitFunction->SetParName(15,"V2"); gFitFunction->SetParameter(15,0.1);
-//   gFitFunction->SetParName(16,"V3"); gFitFunction->SetParameter(16,0.05);
-//   gFitFunction->SetParName(17,"V4"); gFitFunction->SetParameter(17,0.005);
-
-//   // flow parameters
-//   Double_t fNV = 0.;
-//   Double_t fV1 = 0.;
-//   Double_t fV2 = 0.;
-//   Double_t fV3 = 0.;
-//   Double_t fV4 = 0.;
+//____________________________________________________________//
+void drawProjections(const char* lhcPeriod = "LHC10h",
+                    const char* gCentralityEstimator = "V0M",
+                    Int_t gBit = 128,
+                    const char* gEventPlaneEstimator = "VZERO",
+                    Bool_t kProjectInEta = kFALSE,
+                    Int_t binMin = 1,
+                    Int_t binMax = 80,
+                    Int_t gCentrality = 1,
+                    Double_t psiMin = -0.5, 
+                    Double_t psiMax = 3.5,
+                    Double_t vertexZMin = -10., 
+                    Double_t vertexZMax = 10.,
+                    Double_t ptTriggerMin = -1.,
+                    Double_t ptTriggerMax = -1.,
+                    Double_t ptAssociatedMin = -1.,
+                    Double_t ptAssociatedMax = -1.,
+                    Bool_t kUseZYAM = kFALSE) {
+  //Macro that draws the charge dependent correlation functions PROJECTIONS 
+  //for each centrality bin for the different pT of trigger and 
+  //associated particles
+  TGaxis::SetMaxDigits(3);
+
+  //Get the input file
+  /*TString filename = lhcPeriod; 
+  filename += "/Centrality"; filename += gCentralityEstimator;
+  filename += "_Bit"; filename += gBit;
+  filename += "_"; filename += gEventPlaneEstimator;
+  filename +="/PttFrom";
+  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
+  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
+  filename += Form("%.1f",ptAssociatedMax); */
+
+  TString filename = "correlationFunction.Centrality";
+  filename += gCentrality; filename += ".Psi";
+  if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
+  else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
+  else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
+  else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
+  else filename += "All.PttFrom";
+  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
+  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
+  filename += Form("%.1f",ptAssociatedMax); 
+  
+  filename += "_"; 
+  filename += Form("%.1f",psiMin);
+  filename += "-";   
+  filename += Form("%.1f",psiMax);
+  filename += ".root";  
+
+  //Open the file
+  TFile *f = TFile::Open(filename.Data());
+  if((!f)||(!f->IsOpen())) {
+    Printf("The file %s is not found. Aborting...",filename);
+    return listBF;
+  }
+  //f->ls();
+  
+  //Latex
+  TString centralityLatex = "Centrality: ";
+  centralityLatex += centralityArray[gCentrality-1]; 
+  centralityLatex += "%";
+
+  TString psiLatex;
+  if((psiMin == -0.5)&&(psiMax == 0.5))
+    psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}"; 
+  else if((psiMin == 0.5)&&(psiMax == 1.5))
+    psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}"; 
+  else if((psiMin == 1.5)&&(psiMax == 2.5))
+    psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}"; 
+  else 
+    psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}"; 
  
-//   //Fitting the correlation function (first the edges to extract flow)
-//   gHist->Fit("gFitFunction","nm","",1.0,1.6);
-
-//   fNV += gFitFunction->GetParameter(13);
-//   fV1 += gFitFunction->GetParameter(14);
-//   fV2 += gFitFunction->GetParameter(15);
-//   fV3 += gFitFunction->GetParameter(16);
-//   fV4 += gFitFunction->GetParameter(17);
-
-//   gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
-
-//   fNV += gFitFunction->GetParameter(13);
-//   fV1 += gFitFunction->GetParameter(14);
-//   fV2 += gFitFunction->GetParameter(15);
-//   fV3 += gFitFunction->GetParameter(16);
-//   fV4 += gFitFunction->GetParameter(17);
-
-//   // Now fit the whole with fixed flow
-//   gFitFunction->FixParameter(13,fNV/2.);
-//   gFitFunction->FixParameter(14,fV1/2.);
-//   gFitFunction->FixParameter(15,fV2/2.);
-//   gFitFunction->FixParameter(16,fV3/2.);
-//   gFitFunction->FixParameter(17,fV4/2.);
+  TString pttLatex = Form("%.1f",ptTriggerMin);
+  pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
+  pttLatex += " GeV/c";
+
+  TString ptaLatex = Form("%.1f",ptAssociatedMin);
+  ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
+  ptaLatex += " GeV/c";
+
+  TLatex *latexInfo1 = new TLatex();
+  latexInfo1->SetNDC();
+  latexInfo1->SetTextSize(0.045);
+  latexInfo1->SetTextColor(1);
+
+  TString pngName;
+
+  //============================================================//
+  //Get the +- correlation function
+  TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
+  gHistPN->SetStats(kFALSE);
+  gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
+  gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
+  gHistPN->GetXaxis()->CenterTitle();
+  gHistPN->GetXaxis()->SetTitleOffset(1.2);
+  gHistPN->GetYaxis()->CenterTitle();
+  gHistPN->GetYaxis()->SetTitleOffset(1.2);
+  gHistPN->GetZaxis()->SetTitleOffset(1.5);
+  TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
+  cPN->SetFillColor(10); 
+  cPN->SetHighLightColor(10);
+  cPN->SetLeftMargin(0.15);
+
+  //================//
+  TH1D* gHistPNprojection = 0x0;
+  Double_t sum = 0.0;
+  Double_t gError = 0.0;
+  Int_t nCounter = 0;
+  //projection in delta eta
+  if(kProjectInEta) {
+    gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsX(),gHistPN->GetXaxis()->GetXmin(),gHistPN->GetXaxis()->GetXmax());
+    for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
+       sum += gHistPN->GetBinContent(iBinX,iBinY);
+       if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistPNprojection->SetBinContent(iBinX,sum);
+      gHistPNprojection->SetBinError(iBinX,gError);
+    }
+    gHistPNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
+    if(kUseZYAM) 
+      gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
+    else 
+      gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#eta)");
+    gHistPNprojection->GetXaxis()->SetTitle("#Delta#eta");
+  }//projection in delta eta
+  //projection in delta phi
+  else {
+    gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsY(),gHistPN->GetYaxis()->GetXmin(),gHistPN->GetYaxis()->GetXmax());
+    for(Int_t iBinY = 1; iBinY <= gHistPN->GetNbinsY(); iBinY++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
+      //for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
+       sum += gHistPN->GetBinContent(iBinX,iBinY);
+       if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistPNprojection->SetBinContent(iBinY,sum);
+      gHistPNprojection->SetBinError(iBinY,gError);
+    }
+    if(kUseZYAM) 
+      gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
+    else 
+      gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#varphi)");
+    gHistPNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+  }
+
+  //ZYAM
+  if(kUseZYAM) {
+    Double_t reference = gHistPNprojection->GetBinContent(gHistPNprojection->GetMinimumBin());
+    for(Int_t iBinX = 1; iBinX <= gHistPNprojection->GetNbinsX(); iBinX++) 
+      gHistPNprojection->SetBinContent(iBinX,gHistPNprojection->GetBinContent(iBinX) - reference);
+  }
+  
+  gHistPNprojection->GetYaxis()->SetTitleOffset(1.5);
+  gHistPNprojection->SetMarkerStyle(20);
+  gHistPNprojection->SetStats(kFALSE);
+  gHistPNprojection->DrawCopy("E");
+  //=================//
   
-//   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
-//   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
-//   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
-//   gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
-//   gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
-
-//   gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(1,0.3);
-//   gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
-//   gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
-//   gFitFunction->ReleaseParameter(11);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
-//   gFitFunction->ReleaseParameter(12);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
-
-//   gHist->Fit("gFitFunction","nm");
-
-
-//   //Cloning the histogram
-//   TString histoName = gHist->GetName(); histoName += "Fit"; 
-//   TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
-//   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
-//   gHistResidual->SetName("gHistResidual");
-//   gHistResidual->Sumw2();
-
-//   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
-//     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
-//       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
-//     }
-//   }
-//   gHistResidual->Add(gHistFit,-1);
-
-//   //Write to output file
-//   TString newFileName = "correlationFunctionFit";
-//   if(histoName.Contains("PN")) newFileName += "PN";
-//   else if(histoName.Contains("NP")) newFileName += "NP";
-//   else if(histoName.Contains("PP")) newFileName += "PP";
-//   else if(histoName.Contains("NN")) newFileName += "NN";
-//   newFileName += ".Centrality";  
-//   newFileName += gCentrality; newFileName += ".Psi";
-//   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
-//   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
-//   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
-//   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
-//   else newFileName += "All.PttFrom";
-//   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To"; 
-//   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
-//   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
-//   newFileName += Form("%.1f",ptAssociatedMax); 
-//   newFileName += ".root";
-//   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
-//   gHist->Write();
-//   gHistFit->Write();
-//   gHistResidual->Write();
-//   gFitFunction->Write();
-//   newFile->Close();
+  gPad->SetTheta(30); // default is 30
+  gPad->SetPhi(-60); // default is 30
+  gPad->Update();
+
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
+
+  pngName = "Projection.CorrelationFunction.Centrality"; 
+  pngName += centralityArray[gCentrality-1]; 
+  pngName += ".Psi"; 
+  if(kProjectInEta)
+    pngName += ".ETAprojection.";
+  else
+    pngName += ".PHIprojection.";
+  if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+  else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+  else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+  else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+  else pngName += "All.PttFrom";
+  pngName += Form("%.1f",ptTriggerMin); pngName += "To"; 
+  pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+  pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
+  pngName += Form("%.1f",ptAssociatedMax); 
+  pngName += ".PositiveNegative.png";
+  cPN->SaveAs(pngName.Data());
   
+  //============================================================//
+  //Get the -+ correlation function
+  TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
+  gHistNP->SetStats(kFALSE);
+  gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
+  gHistNP->GetXaxis()->SetRangeUser(-1.4,1.4);
+  gHistNP->GetXaxis()->CenterTitle();
+  gHistNP->GetXaxis()->SetTitleOffset(1.2);
+  gHistNP->GetYaxis()->CenterTitle();
+  gHistNP->GetYaxis()->SetTitleOffset(1.2);
+  gHistNP->GetZaxis()->SetTitleOffset(1.5);
+  TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
+  cNP->SetFillColor(10); cNP->SetHighLightColor(10);
+  cNP->SetLeftMargin(0.15);
+
+  //================//
+  TH1D* gHistNPprojection = 0x0;
+  Double_t sum = 0.0;
+  Double_t gError = 0.0;
+  Int_t nCounter = 0;
+  if(kProjectInEta) {
+    gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsX(),gHistNP->GetXaxis()->GetXmin(),gHistNP->GetXaxis()->GetXmax());
+    for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
+      //for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
+       sum += gHistNP->GetBinContent(iBinX,iBinY);
+       if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistNPprojection->SetBinContent(iBinX,sum);
+      gHistNPprojection->SetBinError(iBinX,gError);
+    }
+    gHistNPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
+    if(kUseZYAM) 
+      gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
+    else
+      gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#eta)");
+    gHistNPprojection->GetXaxis()->SetTitle("#Delta#eta");
+  }
+  else {
+    gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsY(),gHistNP->GetYaxis()->GetXmin(),gHistNP->GetYaxis()->GetXmax());
+    for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
+       //for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
+       sum += gHistNP->GetBinContent(iBinX,iBinY);
+       if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistNPprojection->SetBinContent(iBinY,sum);
+      gHistNPprojection->SetBinError(iBinY,gError);
+    }
+    if(kUseZYAM) 
+      gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
+    else
+      gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#varphi)");
+    gHistNPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+  }
+  //ZYAM
+  if(kUseZYAM) {
+    Double_t reference = gHistNPprojection->GetBinContent(gHistNPprojection->GetMinimumBin());
+  for(Int_t iBinX = 1; iBinX <= gHistNPprojection->GetNbinsX(); iBinX++) 
+    gHistNPprojection->SetBinContent(iBinX,gHistNPprojection->GetBinContent(iBinX) - reference);
+  }
+
+  gHistNPprojection->GetYaxis()->SetTitleOffset(1.5);
+  gHistNPprojection->SetMarkerStyle(20);
+  gHistNPprojection->SetStats(kFALSE);
+  gHistNPprojection->DrawCopy("E");
+  //================//
+
+  gPad->SetTheta(30); // default is 30
+  gPad->SetPhi(-60); // default is 30
+  gPad->Update();
+
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
+
+  pngName = "Projection.CorrelationFunction.Centrality"; 
+  pngName += centralityArray[gCentrality-1]; 
+  pngName += ".Psi"; 
+  if(kProjectInEta)
+    pngName += ".ETAprojection.";
+  else
+    pngName += ".PHIprojection.";
+  if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+  else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+  else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+  else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+  else pngName += "All.PttFrom";
+  pngName += Form("%.1f",ptTriggerMin); pngName += "To"; 
+  pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+  pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
+  pngName += Form("%.1f",ptAssociatedMax); 
+  pngName += ".NegativePositive.png";
+  cNP->SaveAs(pngName.Data());
+
+  //============================================================//
+  //Get the ++ correlation function
+  TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
+  gHistPP->SetStats(kFALSE);
+  gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
+  gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
+  gHistPP->GetXaxis()->CenterTitle();
+  gHistPP->GetXaxis()->SetTitleOffset(1.2);
+  gHistPP->GetYaxis()->CenterTitle();
+  gHistPP->GetYaxis()->SetTitleOffset(1.2);
+  gHistPP->GetZaxis()->SetTitleOffset(1.5);
+  TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
+  cPP->SetFillColor(10); cPP->SetHighLightColor(10);
+  cPP->SetLeftMargin(0.15);
+  //=================//
+  TH1D* gHistPPprojection = 0x0;
+  Double_t sum = 0.0;
+  Double_t gError = 0.0;
+  Int_t nCounter = 0;
+  if(kProjectInEta) {
+    gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsX(),gHistPP->GetXaxis()->GetXmin(),gHistPP->GetXaxis()->GetXmax());
+    for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
+       //for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
+       sum += gHistPP->GetBinContent(iBinX,iBinY);
+       if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistPPprojection->SetBinContent(iBinX,sum);
+      gHistPPprojection->SetBinError(iBinX,gError);
+    }
+    gHistPPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
+    if(kUseZYAM) 
+      gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
+    else
+      gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#eta)");
+    gHistPPprojection->GetXaxis()->SetTitle("#Delta#eta");
+  }
+  else {
+    gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsY(),gHistPP->GetYaxis()->GetXmin(),gHistPP->GetYaxis()->GetXmax());
+    for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
+       //for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
+       sum += gHistPP->GetBinContent(iBinX,iBinY);
+       if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistPPprojection->SetBinContent(iBinY,sum);
+      gHistPPprojection->SetBinError(iBinY,gError);
+    }
+    if(kUseZYAM) 
+      gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
+    else
+      gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#varphi)");
+    gHistPPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+  }
+  //ZYAM
+  if(kUseZYAM) {
+    Double_t reference = gHistPPprojection->GetBinContent(gHistPPprojection->GetMinimumBin());
+    for(Int_t iBinX = 1; iBinX <= gHistPPprojection->GetNbinsX(); iBinX++) 
+      gHistPPprojection->SetBinContent(iBinX,gHistPPprojection->GetBinContent(iBinX) - reference);
+  }
+
+  gHistPPprojection->GetYaxis()->SetTitleOffset(1.5);
+  gHistPPprojection->SetMarkerStyle(20);
+  gHistPPprojection->SetStats(kFALSE);
+  gHistPPprojection->DrawCopy("E");
+  //================//
+
+  gPad->SetTheta(30); // default is 30
+  gPad->SetPhi(-60); // default is 30
+  gPad->Update();
+
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
+
+  pngName = "Projection.CorrelationFunction.Centrality"; 
+  pngName += centralityArray[gCentrality-1]; 
+  pngName += ".Psi"; 
+  if(kProjectInEta) 
+    pngName += ".ETAprojection.";
+  else
+    pngName += ".PHIprojection.";
+  if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+  else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+  else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+  else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+  else pngName += "All.PttFrom";
+  pngName += Form("%.1f",ptTriggerMin); pngName += "To"; 
+  pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+  pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
+  pngName += Form("%.1f",ptAssociatedMax); 
+  pngName += ".PositivePositive.png";
+  cPP->SaveAs(pngName.Data());
+
+  //============================================================//
+  //Get the -- correlation function
+  TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
+  gHistNN->SetStats(kFALSE);
+  gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
+  gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
+  gHistNN->GetXaxis()->CenterTitle();
+  gHistNN->GetXaxis()->SetTitleOffset(1.2);
+  gHistNN->GetYaxis()->CenterTitle();
+  gHistNN->GetYaxis()->SetTitleOffset(1.2);
+  gHistNN->GetZaxis()->SetTitleOffset(1.5);
+  TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
+  cNN->SetFillColor(10); cNN->SetHighLightColor(10);
+  cNN->SetLeftMargin(0.15);
+
+  //=================//
+  TH1D* gHistNNprojection = 0x0;
+  Double_t sum = 0.0;
+  Double_t gError = 0.0;
+  Int_t nCounter = 0;
+  if(kProjectInEta) {
+    gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsX(),gHistNN->GetXaxis()->GetXmin(),gHistNN->GetXaxis()->GetXmax());
+    for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
+       //for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
+       sum += gHistNN->GetBinContent(iBinX,iBinY);
+       if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistNNprojection->SetBinContent(iBinX,sum);
+      gHistNNprojection->SetBinError(iBinX,gError);
+    }
+    gHistNNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
+    if(kUseZYAM) 
+      gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
+    else
+      gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#eta)");
+    gHistNNprojection->GetXaxis()->SetTitle("#Delta#eta");
+  }
+  else {
+    gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsY(),gHistNN->GetYaxis()->GetXmin(),gHistNN->GetYaxis()->GetXmax());
+    for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
+      sum = 0.; gError = 0.0; nCounter = 0;
+      for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
+       //for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
+       sum += gHistNN->GetBinContent(iBinX,iBinY);
+       if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
+        Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
+       gError += exy*exy;
+      }
+      if(nCounter != 0) {
+       sum /= nCounter;
+       gError = TMath::Sqrt(gError)/nCounter;
+      }
+      gHistNNprojection->SetBinContent(iBinY,sum);
+      gHistNNprojection->SetBinError(iBinY,gError);
+    }
+    if(kUseZYAM) 
+      gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
+    else
+      gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#varphi)");
+    gHistNNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+  }
+  //ZYAM
+  if(kUseZYAM) {
+    Double_t reference = gHistNNprojection->GetBinContent(gHistNNprojection->GetMinimumBin());
+    for(Int_t iBinX = 1; iBinX <= gHistNNprojection->GetNbinsX(); iBinX++) 
+      gHistNNprojection->SetBinContent(iBinX,gHistNNprojection->GetBinContent(iBinX) - reference); 
+  }
+
+  gHistNNprojection->GetYaxis()->SetTitleOffset(1.5);
+  gHistNNprojection->SetMarkerStyle(20);
+  gHistNNprojection->SetStats(kFALSE);
+  gHistNNprojection->DrawCopy("E");
+  //=================//
+
+  gPad->SetTheta(30); // default is 30
+  gPad->SetPhi(-60); // default is 30
+  gPad->Update();
+
+  latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
+  latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
+
+  pngName = "Projection.CorrelationFunction.Centrality"; 
+  pngName += centralityArray[gCentrality-1]; 
+  pngName += ".Psi"; 
+  if(kProjectInEta) 
+    pngName += ".ETAprojection.";
+  else
+    pngName += ".PHIprojection.";
+  if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+  else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+  else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+  else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+  else pngName += "All.PttFrom";
+  pngName += Form("%.1f",ptTriggerMin); pngName += "To"; 
+  pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+  pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
+  pngName += Form("%.1f",ptAssociatedMax); 
+  pngName += ".NegativeNegative.png";
+  cNN->SaveAs(pngName.Data());
 
-// }
+  TFile *fProjection = TFile::Open("ProjectionCorrelationFunction.root",
+                                         "recreate");
+  
+  gHistNPprojection->Write();
+  gHistPNprojection->Write();
+  gHistNNprojection->Write();
+  gHistPPprojection->Write();
+  fProjection->Close();
+}
 
 //____________________________________________________________//
 void fitCorrelationFunctions(Int_t gCentrality = 1,
                             Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+                            Double_t vertexZMin = -10.,Double_t vertexZMax = 10.,
                             Double_t ptTriggerMin = -1.,
                             Double_t ptTriggerMax = -1.,
                             Double_t ptAssociatedMin = -1.,
@@ -1864,6 +2254,12 @@ void fitCorrelationFunctions(Int_t gCentrality = 1,
   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
   newFileName += Form("%.1f",ptAssociatedMax); 
+
+  newFileName += "_"; 
+  newFileName += Form("%.1f",psiMin);
+  newFileName += "-";   
+  newFileName += Form("%.1f",psiMax); 
+
   newFileName += ".root";
   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
   gHist->Write();
@@ -1875,105 +2271,8 @@ void fitCorrelationFunctions(Int_t gCentrality = 1,
 
 }
 
-
-
-// //____________________________________________________________//
-// void fitCorrelationFunctions(Int_t gCentrality = 1,
-//                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
-//                          Double_t ptTriggerMin = -1.,
-//                          Double_t ptTriggerMax = -1.,
-//                          Double_t ptAssociatedMin = -1.,
-//                          Double_t ptAssociatedMax = -1.,
-//                          TH2D *gHist) {
-
-//   cout<<"FITTING FUNCTION (HOUSTON)"<<endl;
-
-//   // Fit Function
-//    //x axis = delta_eta
-//    //y axis = delta_phi
-//                                                                                                                                                                                         //  [9]*exp(-1*pow(((x/[10])^2 + (y/[10])^2),0.5)) Hyper expo
-//   TF2 *fit1 = new TF2("fit1","[0] + [1]*cos(y) + [2]*cos(2*y) + [3]*cos(3*y) + [4]*cos(4*y) +[5]*cos(5*y)+ [6]*exp(-0.5*pow(((x/[7])^2 + (y/[8])^2),[11])) + [6]*exp(-0.5*pow(((x/[7])^2 + ((y-6.283)/[8])^2),[11]))+ [9]*exp(-1*((x/[10])^2 + (y/[10])^2)) ",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
-  
-   
-//   Double_t Parameters[] = {0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.3,0.3,1.0,0.1,0.1};
-       
-//   fit1->SetParameters(Parameters);  // input pars from macro arguments
-       
-//   fit1->SetParName(0,"offset");
-//   fit1->SetParName(1,"v1");
-//   fit1->SetParName(2,"v2");
-//   fit1->SetParName(3,"v3");
-//   fit1->SetParName(4,"v4");
-//   fit1->SetParName(5,"v5");
-//   fit1->SetParName(6,"Ridgeamp");
-//   fit1->SetParName(7,"Ridgesigx");
-//   fit1->SetParName(8,"Ridgesigy");
-//   fit1->SetParName(9,"Expoamp");
-//   fit1->SetParName(10,"Exposig");
-//   fit1->SetParName(11,"Gausspara");
-    
-
-//   //Fit Parameter Ranges
-//   fit1->SetParLimits(0,-2.0,2.0);   //offset
-//   fit1->SetParLimits(1,-1.0,0.1);   //v1
-//   fit1->SetParLimits(2,-1.6,0.9);   //v2
-//   fit1->SetParLimits(3,0.0,0.5);    //v3
-//   fit1->SetParLimits(4,0.0,0.9);    //v4
-//   fit1->SetParLimits(5,0.0,0.9);    //v5
-//   fit1->SetParLimits(6,0.0,1.5);    //Ridgeamp
-//   fit1->SetParLimits(7,0.1,3.0);    //Ridgesigx
-//   fit1->SetParLimits(8,0.1,2.0);    //Ridgesigy
-//   fit1->SetParLimits(9,0.0,6.0);      //Expoamp
-//   fit1->SetParLimits(10,0.05,0.5); //Exposig
-//   fit1->SetParLimits(11,0.0,2.0);    //Gausspara
-
-//   //Fitting the correlation function
-//   gHist->Fit("fit1","nm");
-
-//   //Cloning the histogram
-//   TString histoName = gHist->GetName(); histoName += "Fit"; 
-//   TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
-//   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
-//   gHistResidual->SetName("gHistResidual");
-//   gHistResidual->Sumw2();
-
-//   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
-//     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
-//       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,fit1->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
-//     }
-//   }
-//   gHistResidual->Add(gHistFit,-1);
-
-//   //Write to output file
-//   TString newFileName = "correlationFunctionFit";
-//   if(histoName.Contains("PN")) newFileName += "PN";
-//   else if(histoName.Contains("NP")) newFileName += "NP";
-//   else if(histoName.Contains("PP")) newFileName += "PP";
-//   else if(histoName.Contains("NN")) newFileName += "NN";
-//   newFileName += ".Centrality";  
-//   newFileName += gCentrality; newFileName += ".Psi";
-//   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
-//   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
-//   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
-//   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
-//   else newFileName += "All.PttFrom";
-//   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To"; 
-//   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
-//   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To"; 
-//   newFileName += Form("%.1f",ptAssociatedMax); 
-//   newFileName += ".root";
-//   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
-//   gHist->Write();
-//   gHistFit->Write();
-//   gHistResidual->Write();
-//   fit1->Write();
-//   newFile->Close();
-  
-
-// }
-
-TH2D *convolute2D(TH2D* h1, TH2D* h2, TString hname){
+//____________________________________________________________//
+TH2D *convolute2D(TH2D* h1, TH2D* h2, TString hname) {
   //
   // this function does the convolution of 2 eta or phi "efficiencies" in a deta or dphi "efficiency"
   // and returns a new histogram which is normalized to the number of bin combinations