Getting the bf from correlation functions

diff --git a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
index a4e856c084dc6771b0bc38c9b8d39f10ba6e60bd..a42dbd32e4c901132e0d549e075bf406b6fb7144 100644 (file)
--- a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
+++ b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
@@ -908,7 +908,8 @@ void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
 }
 
 //____________________________________________________________//
-void drawProjections(Bool_t kProjectInEta = kFALSE,
+void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
+                    Bool_t kProjectInEta = kFALSE,
                     Int_t binMin = 1,
                     Int_t binMax = 80,
                     Int_t gCentrality = 1,
@@ -978,13 +979,16 @@ void drawProjections(Bool_t kProjectInEta = kFALSE,
   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;
+  TFile *f = 0x0;
+  if(!gHistBalanceFunction2D) {
+    TFile::Open(filename.Data());
+    if((!f)||(!f->IsOpen())) {
+      Printf("The file %s is not found. Aborting...",filename);
+      return listBF;
+    }
+    //f->ls();
   }
-  //f->ls();
-  
+
   //Latex
   TString centralityLatex = "Centrality: ";
   centralityLatex += centralityArray[gCentrality-1]; 
@@ -1016,9 +1020,16 @@ void drawProjections(Bool_t kProjectInEta = kFALSE,
 
   //============================================================//
   //Get subtracted and mixed balance function
-
-  TH2D *gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
-  TH2D *gHistBalanceFunctionMixed2D      = (TH2D*)f->Get("gHistBalanceFunctionMixed");
+  TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
+  TH2D *gHistBalanceFunctionMixed2D      = 0x0;
+  if(!gHistBalanceFunction2D) {
+    gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
+    gHistBalanceFunctionMixed2D      = (TH2D*)f->Get("gHistBalanceFunctionMixed");
+  }
+  else {
+    gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
+    gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
+  }
 
   TH1D *gHistBalanceFunctionSubtracted = NULL;
   TH1D *gHistBalanceFunctionMixed      = NULL;
@@ -1098,35 +1109,49 @@ void drawProjections(Bool_t kProjectInEta = kFALSE,
 
 
     // store in txt files
-
     TString meanFileName = filename;
-    if(kProjectInEta) meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
-    else              meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
-    ofstream fileMean(meanFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      meanFileName= "deltaEtaProjection_Mean.txt";
+    //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
+    else              
+      meanFileName = "deltaPhiProjection_Mean.txt";
+      //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
+    ofstream fileMean(meanFileName.Data(),ios::app);
     fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
     fileMean.close();
 
     TString rmsFileName = filename;
-    if(kProjectInEta) rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
-    else              rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
-    ofstream fileRms(rmsFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      rmsFileName = "deltaEtaProjection_Rms.txt";
+      //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
+    else              
+      rmsFileName = "deltaPhiProjection_Rms.txt");
+      //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
+    ofstream fileRms(rmsFileName.Data(),ios::app);
     fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
     fileRms.close();
 
     TString skewnessFileName = filename;
-    if(kProjectInEta) skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
-    else              skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
-    ofstream fileSkewness(skewnessFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      skewnessFileName = "deltaEtaProjection_Skewness.txt";
+      //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
+    else              
+      skewnessFileName = "deltaPhiProjection_Skewness.txt";
+    //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
+    ofstream fileSkewness(skewnessFileName.Data(),ios::app);
     fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
     fileSkewness.close();
 
     TString kurtosisFileName = filename;
-    if(kProjectInEta) kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
-    else              kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
-    ofstream fileKurtosis(kurtosisFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
+      //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
+    else
+      kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";              
+      //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
+    ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
     fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
     fileKurtosis.close();
-   
   }
   // calculate the moments by hand
   else{
@@ -1168,30 +1193,46 @@ void drawProjections(Bool_t kProjectInEta = kFALSE,
 
     // store in txt files
     TString meanFileName = filename;
-    if(kProjectInEta) meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
-    else              meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
-    ofstream fileMean(meanFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      meanFileName = "deltaEtaProjection_Mean.txt";
+      //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
+    else              
+      meanFileName = "deltaPhiProjection_Mean.txt";
+      //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
+    ofstream fileMean(meanFileName.Data(),ios::app);
     fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
     fileMean.close();
 
     TString rmsFileName = filename;
-    if(kProjectInEta) rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
-    else              rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
-    ofstream fileRms(rmsFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      rmsFileName = "deltaEtaProjection_Rms.txt";
+      //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
+    else              
+      rmsFileName = "deltaPhiProjection_Rms.txt";
+//rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
+    ofstream fileRms(rmsFileName.Data(),ios::app);
     fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
     fileRms.close();
 
     TString skewnessFileName = filename;
-    if(kProjectInEta) skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
-    else              skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
-    ofstream fileSkewness(skewnessFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      skewnessFileName = "deltaEtaProjection_Skewness.txt";
+      //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
+    else              
+      skewnessFileName = "deltaPhiProjection_Skewness.txt";
+      //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
+    ofstream fileSkewness(skewnessFileName.Data(),ios::app);
     fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
     fileSkewness.close();
 
     TString kurtosisFileName = filename;
-    if(kProjectInEta) kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
-    else              kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
-    ofstream fileKurtosis(kurtosisFileName.Data(),ios::out);
+    if(kProjectInEta) 
+      kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
+      //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
+    else              
+      kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
+      //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
+    ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
     fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
     fileKurtosis.close();
   }
@@ -1227,3 +1268,143 @@ void drawProjections(Bool_t kProjectInEta = kFALSE,
   gHistBalanceFunctionMixed->Write();
   fProjection->Close();
 }
+
+//____________________________________________________________//
+void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
+                                        Int_t gTrainNumber = 64,
+                                        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.) {
+  //Macro that draws the BF distributions for each centrality bin
+  //for reaction plane dependent analysis
+  //Author: Panos.Christakoglou@nikhef.nl
+  TGaxis::SetMaxDigits(3);
+  gStyle->SetPalette(55,0);
+
+  //Get the input file
+  TString filename = lhcPeriod; 
+  //filename += "/Train"; filename += gTrainNumber;
+  filename +="/PttFrom";
+  filename += Form("%.1f",ptTriggerMin); filename += "To"; 
+  filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+  filename += Form("%.1f",ptAssociatedMin); filename += "To"; 
+  filename += Form("%.1f",ptAssociatedMax); 
+  filename += "/correlationFunction.";
+  filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
+  filename += ".PsiAll.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();
+
+  TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
+  if(!gHistPN) return;
+  TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
+  if(!gHistNP) return;
+  TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
+  if(!gHistPP) return;
+  TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
+  if(!gHistNN) return;
+
+  gHistPN->Sumw2();
+  gHistPP->Sumw2();
+  gHistPN->Add(gHistPP,-1);
+  gHistNP->Sumw2();
+  gHistNN->Sumw2();
+  gHistNP->Add(gHistNN,-1);
+  gHistPN->Add(gHistNP);
+  gHistPN->Scale(0.5);
+  TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
+  gHistBalanceFunction2D->SetStats(kFALSE);
+  gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
+  gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
+  gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
+
+  //Draw the results
+  TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
+  c0->SetFillColor(10); c0->SetHighLightColor(10);
+  c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
+  gHistBalanceFunction2D->SetTitle("");
+  gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
+  gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
+  gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
+  gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
+  gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4); 
+  gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
+  gHistBalanceFunction2D->DrawCopy("lego2");
+  gPad->SetTheta(30); // default is 30
+  gPad->SetPhi(-60); // default is 30
+  gPad->Update();  
+ 
+  TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
+
+  TString pttLatex = Form("%.1f",ptTriggerMin);
+  pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
+  pttLatex += " GeV/c";
+
+  TString ptaLatex = Form("%.1f",ptAssociatedMin);
+  ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
+  ptaLatex += " GeV/c";
+
+  TLatex *latexInfo1 = new TLatex();
+  latexInfo1->SetNDC();
+  latexInfo1->SetTextSize(0.045);
+  latexInfo1->SetTextColor(1);
+  latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
+  latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
+  latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
+
+  TString pngName = "BalanceFunction2D."; 
+  pngName += Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
+  pngName += ".PttFrom";  
+  pngName += Form("%.1f",ptTriggerMin); pngName += "To"; 
+  pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+  pngName += Form("%.1f",ptAssociatedMin); pngName += "To"; 
+  pngName += Form("%.1f",ptAssociatedMax); 
+  pngName += ".png";
+
+  c0->SaveAs(pngName.Data());
+
+  drawProjections(gHistBalanceFunction2D,
+                 kTRUE,
+                 1,36,
+                 gCentrality,
+                 psiMin,psiMax,
+                 ptTriggerMin,ptTriggerMax,
+                 ptAssociatedMin,ptAssociatedMax,
+                 kTRUE,
+                 "Multiplicity",
+                 kFALSE);
+
+  drawProjections(gHistBalanceFunction2D,
+                 kFALSE,
+                 1,80,
+                 gCentrality,
+                 psiMin,psiMax,
+                 ptTriggerMin,ptTriggerMax,
+                 ptAssociatedMin,ptAssociatedMax,
+                 kTRUE,
+                 "Multiplicity",
+                 kFALSE);
+}