Adding FWHM to determine BF width (Alis Rodriguez Manso <alisrm@nikhef.nl>)

diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
index 011c640fbeb5d79c06d12f347b02bcc99dcea236..f93257604e368e51c3ccd22b6247c0fd1323175a 100644 (file)
--- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
+++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
@@ -31,6 +31,8 @@
 #include <TObjArray.h>
 #include <TGraphErrors.h>
 #include <TString.h>
+#include <TSpline.h>
+#include <TRandom3.h>
 
 #include "AliVParticle.h"
 #include "AliMCParticle.h"
@@ -2431,7 +2433,6 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
 }
 
 //____________________________________________________________________//
-
 Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
                                           Double_t &mean, Double_t &meanError,
                                           Double_t &sigma, Double_t &sigmaError,
@@ -2654,3 +2655,223 @@ Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag,
   AliFatal(Form("Tag %s not found in %s", tag, configuration));
   return 0;
 }
+
+//____________________________________________________________________//
+Double_t AliBalancePsi::GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist,
+                               Double_t &fwhm_spline, Double_t &fwhmError) {
+ 
+  // helper method to calculate the fwhm and errors of a TH1D 
+  if(gHist){
+    Int_t repeat = 10000;
+    
+    if (gDeltaEtaPhi == 1){ 
+      fwhm_spline  = 0.;
+      fwhmError    = 0.;
+      gHist->Smooth(4); 
+
+      Double_t xmin = gHist->GetXaxis()->GetXmin();
+      Double_t xmax = gHist->GetXaxis()->GetXmax();
+      
+      //+++++++FWHM TSpline+++++++++++++++++++++++++++//
+      TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); 
+      spline3->SetLineColor(kGreen+2); 
+      spline3->SetLineWidth(2);
+      //spline3->Draw("SAME");
+      
+      Int_t Nbin = gHist->GetNbinsX();
+      Int_t bin_max_hist_y = gHist->GetMaximumBin();
+      Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
+      Double_t max_spline = spline3->Eval(bins_center_x);      
+      Double_t y_spline = -999;
+      Double_t y_spline_r = -999;
+      Double_t x_spline_l = -999;
+      Double_t x_spline_r = -999;
+      
+      for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
+       y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
+       if (y_spline > max_spline/2){
+         x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
+         break;
+       }
+      } 
+      for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
+       y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
+       if (y_spline_r < max_spline/2){
+         x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
+         break;
+       }
+      }   
+      fwhm_spline = x_spline_r - x_spline_l;
+      //+++++++FWHM TSpline++++++++++++++++++++++++//
+    
+      //+++++++Error Calculation SPLINE++++++++++++//
+      Double_t dmeans,dsigmas;
+      TSpline3 *spline1;  
+      Int_t bin_max_hist_y1;
+      Double_t bins_center_x1;
+      Double_t max_spline1;
+      Double_t y_spline1 = - 999.; 
+      Double_t y_spline_r1 = - 999.; 
+      Double_t x_spline_l1 = -999.; 
+      Double_t x_spline_r1 = -999.;
+      Double_t fwhm_spline1 = 0.;
+      Double_t fwhm_T = 0.;  
+ 
+      TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
+      TRandom3 *rgauss = new TRandom3(0);      
+      TH1F *hists = new TH1F("hists","hists",1000,1.,3);
+      
+      for (Int_t f=0; f<repeat; f++){ //100
+       for (Int_t h=0; h<Nbin+1; h++){   
+         dmeans = gHist->GetBinContent(h);
+         dsigmas = gHist->GetBinError(h);
+         Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
+         hEmpty2->SetBinContent(h,rgauss_point);
+         
+         //++++++++++++//  
+         hEmpty2->Smooth(4);   
+         //++++++++++++//
+       }           
+       spline1 = new TSpline3(hEmpty2,"b2e2",0,0); 
+       spline1->SetLineColor(kMagenta+1); 
+       spline1->SetLineWidth(1);
+       //spline1->Draw("SAME");  
+       
+       bin_max_hist_y1 = hEmpty2->GetMaximumBin();
+       bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
+       max_spline1 = spline1->Eval(bins_center_x1);
+       
+       for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
+         y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
+         if (y_spline1 > max_spline1/2){
+           x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
+           break;
+         }
+       } 
+       for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
+         y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
+         if (y_spline_r1 < max_spline1/2){
+           x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
+           break;
+         }
+       }     
+       
+       fwhm_spline1 = x_spline_r1 - x_spline_l1;
+       hists->Fill(fwhm_spline1);   
+       fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1); 
+      }
+     
+      fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
+      //+++++++Error Calculation SPLINE+++++++++++//
+    }
+    else{  
+      fwhm_spline  = 0.;
+      fwhmError    = 0.;
+      gHist->Smooth(4);
+      
+      Double_t xmin = gHist->GetXaxis()->GetXmin();
+      Double_t xmax = gHist->GetXaxis()->GetXmax();
+      
+      //+++++++FWHM TSpline+++++++++++++++++++++++++++//
+      TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0); 
+      spline3->SetLineColor(kGreen+2); 
+      spline3->SetLineWidth(2);
+      //spline3->Draw("SAME");  
+      
+      Int_t Nbin = gHist->GetNbinsX();
+      Int_t bin_max_hist_y = gHist->GetMaximumBin();
+      Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
+      Double_t max_spline = spline3->Eval(bins_center_x);
+      
+      Int_t bin_min_hist_y = gHist->GetMinimumBin();
+      Double_t bins_center_xmin = gHist->GetBinCenter(bin_min_hist_y);
+      Double_t min_spline = spline3->Eval(bins_center_xmin);
+      
+      Double_t y_spline = -999.;
+      Double_t y_spline_r = -999.;
+      Double_t x_spline_l = -999.;
+      Double_t x_spline_r = -999.;
+      
+      for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
+       y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
+       if (y_spline > (max_spline+min_spline)/2){
+         x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
+         break;
+       }
+      } 
+      for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
+       y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
+       if (y_spline_r < (max_spline+min_spline)/2){
+         x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
+         break;
+       }
+      }   
+      fwhm_spline = x_spline_r - x_spline_l;
+      //+++++++FWHM TSpline++++++++++++++++++++++++//
+      
+      //+++++++Error Calculation SPLINE++++++++++++//
+      Double_t dmeans,dsigmas;
+      TSpline3 *spline1;  
+      Int_t bin_max_hist_y1;
+      Double_t bins_center_x1;
+      Double_t max_spline1;
+      Double_t y_spline1 = -999.;
+      Double_t y_spline_r1 = -999.;
+      Double_t x_spline_l1 = -999.;
+      Double_t x_spline_r1 = -999.;
+      Double_t fwhm_spline1 = 0.;
+      Double_t fwhm_T = 0.;   
+
+      TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
+      TRandom3 *rgauss = new TRandom3(0);      
+      Int_t bin_min_hist_y1; 
+      Double_t bins_center_xmin1,min_spline1;
+      
+      for (Int_t f=0; f<repeat; f++){ //100
+       for (Int_t h=0; h<Nbin+1; h++){   
+         dmeans = gHist->GetBinContent(h);
+         dsigmas = gHist->GetBinError(h);
+         Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
+         hEmpty2->SetBinContent(h,rgauss_point);
+         
+         //++++++++++++//  
+         hEmpty2->Smooth(4);   
+         //++++++++++++//
+       }           
+       spline1 = new TSpline3(hEmpty2,"b2e2",0,0); 
+       spline1->SetLineColor(kMagenta+1); 
+       spline1->SetLineWidth(1);
+       //spline1->Draw("SAME");  
+       
+       bin_max_hist_y1 = hEmpty2->GetMaximumBin();
+       bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
+       max_spline1 = spline1->Eval(bins_center_x1);
+       
+       bin_min_hist_y1 = hEmpty2->GetMinimumBin();
+       bins_center_xmin1 = hEmpty2->GetBinCenter(bin_min_hist_y1);
+       min_spline1 = spline1->Eval(bins_center_xmin1);
+       
+       for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
+         y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
+         if (y_spline1 > (max_spline1+min_spline1)/2){
+           x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
+           break;
+         }
+       } 
+       for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
+         y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
+         if (y_spline_r1 < (max_spline1+min_spline1)/2){
+           x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
+           break;
+         }
+       }     
+       
+       fwhm_spline1 = x_spline_r1 - x_spline_l1; 
+       fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);    
+      }           
+      fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
+    }
+  }
+  return fwhm_spline;
+  return fwhmError;
+}

diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
index dcd7683f0d18b18d97b7e1d7099571a4370d3caa..27f1448513e8405e935e5d80680b85834620f339 100644 (file)
--- a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
+++ b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
@@ -193,6 +193,11 @@ class AliBalancePsi : public TObject {
                              Double_t &sigma, Double_t &sigmaError,
                              Double_t &skewness, Double_t &skewnessError,
                              Double_t &kurtosis, Double_t &kurtosisError);
+
+  //++++++++++++++++++//
+  Double_t GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist,
+                  Double_t &fwhm_spline, Double_t &fwhmError);
+  //++++++++++++++++++//
   
   TH2D *GetQAHistHBTbefore() {return fHistHBTbefore;}
   TH2D *GetQAHistHBTafter() {return fHistHBTafter;}

diff --git a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
index d8a372940447f358fb79ed20141f948973ffd63b..5613ff5449cb744b22aecc9307d3c22bb06660bd 100644 (file)
--- a/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
+++ b/PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
@@ -1,5 +1,5 @@
 const Int_t numberOfCentralityBins = 12;
-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"};
+TString centralityArray[numberOfCentralityBins] = {"0-80","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;
@@ -933,7 +933,8 @@ void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
                     TString eventClass = "Centrality",
                     Bool_t bRootMoments = kFALSE,
                     Bool_t kUseZYAM = kFALSE,
-                    Bool_t bReduceRangeForMoments = kFALSE) {
+                    Bool_t bReduceRangeForMoments = kFALSE,
+                    Bool_t bFWHM = kFALSE) {
   //Macro that draws the balance functions PROJECTIONS 
   //for each centrality bin for the different pT of trigger and 
   //associated particles
@@ -1292,12 +1293,37 @@ void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
   fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
   fileWeightedMean.close();
 
-
   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
   c2->SetFillColor(10); 
   c2->SetHighLightColor(10);
   c2->SetLeftMargin(0.15);
   gHistBalanceFunctionSubtracted->DrawCopy("E");
+
+  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
+  Double_t fwhm_spline = 0.;
+  Double_t fwhmError = 0.;
+  if (bFWHM){ 
+    AliBalancePsi *bHelper_1 = new AliBalancePsi;
+    bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError);
+    Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError);
+    
+    //store and in txt files FWHM
+    TString sigmaFileName = filename;
+    if(kProjectInEta) {
+      sigmaFileName = "deltaEtaProjection_FWHM.txt"; 
+    }
+    else{
+      sigmaFileName = "deltaPhiProjection_FWHM.txt"; 
+    }    
+    ofstream fileSigma(sigmaFileName.Data(),ios::app);
+    fileSigma << " " << fwhm_spline << " " <<fwhmError<<endl;
+    fileSigma.close();
+  
+    gHistBalanceFunctionSubtracted->SetLineColor(2);
+    gHistBalanceFunctionSubtracted->SetMarkerColor(2);
+    gHistBalanceFunctionSubtracted->Draw("SAME");
+  }
+  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
  
   TLatex *latex = new TLatex();
   latex->SetNDC();
@@ -1340,8 +1366,8 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
                                         TString eventClass = "Multiplicity",
                                         Bool_t bRootMoments = kFALSE,
                                         Bool_t bFullPhiForEtaProjection = kFALSE,
-                                        Bool_t bReduceRangeForMoments = kFALSE
-) {
+                                        Bool_t bReduceRangeForMoments = kFALSE,
+                                        Bool_t bFWHM = kFALSE) {
   //Macro that draws the BF distributions for each centrality bin
   //for reaction plane dependent analysis
   //Author: Panos.Christakoglou@nikhef.nl
@@ -1479,10 +1505,11 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
                   eventClass.Data(),
                   bRootMoments,
                   kFALSE,
-                  bReduceRangeForMoments);
+                  bReduceRangeForMoments,
+                  bFWHM);
  }
  else{
-  drawProjections(gHistBalanceFunction2D,
+   drawProjections(gHistBalanceFunction2D,
                   kTRUE,
                   1,36,
                   gCentrality,
@@ -1491,9 +1518,10 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
                   ptAssociatedMin,ptAssociatedMax,
                   kTRUE,
                   eventClass.Data(),
-                 bRootMoments,
-                 kFALSE,
-                 bReduceRangeForMoments);
+                  bRootMoments,
+                  kFALSE,
+                  bReduceRangeForMoments,
+                  bFWHM);
  }
 
   drawProjections(gHistBalanceFunction2D,
@@ -1507,7 +1535,8 @@ void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
                  eventClass.Data(),
                  bRootMoments,
                  kFALSE,
-                 bReduceRangeForMoments);
+                 bReduceRangeForMoments,
+                 bFWHM);
 
   TString outFileName = filename;
   outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
@@ -1643,17 +1672,18 @@ sFileName[iDir] += momDirectory;
  hBFOut->Scale(1./entriesOut);
 
  drawProjections(hBFOut,
-                 kTRUE,
-                 1,36,
-                 gCentrality,
-                 psiMin,psiMax,
-                 ptTriggerMin,ptTriggerMax,
-                 ptAssociatedMin,ptAssociatedMax,
-                 kTRUE,
-                 eventClass,
-                 kTRUE,
-                 kFALSE,
-                 bReduceRangeForMoments);
+                kTRUE,
+                1,36,
+                gCentrality,
+                psiMin,psiMax,
+                ptTriggerMin,ptTriggerMax,
+                ptAssociatedMin,ptAssociatedMax,
+                kTRUE,
+                eventClass,
+                kTRUE,
+                kFALSE,
+                bReduceRangeForMoments,
+                bFWHM);
 
   drawProjections(hBFOut,
                  kFALSE,
@@ -1666,7 +1696,8 @@ sFileName[iDir] += momDirectory;
                  eventClass.Data(),
                  kTRUE,
                  kFALSE,
-                 bReduceRangeForMoments);
+                 bReduceRangeForMoments,
+                 bFWHM);
 
  // output
  TString outFileName = "balanceFunction2D.";